This document provides a step-by-step tutorial on how to perform a remoteness analysis for a user-defined area using R, QGIS and Google Earth Engine. Remoteness indicates the distance to roads, but takes into account the cost to walk through terrain. It can serve as anthropogenic factor in wildlife studies, where it has been shown to be an important driver of species occurrence. MORE ON APPLICATIONS?…
Input data for the remoteness analysis need to be prepared in R and
QGIS. The remoteness analysis needs to be performed in Google Earth
Engine.
Using the province of Kâmpóng Thum in Cambodia as an example area,
this tutorial presents in detail all steps required for the
analysis.
Figure: Final remoteness layer calculated from osm drivable
roads and applied water mask.
Useful information when performing remoteness analysis:
# Set path to working directory
wd <- "D:/Dateien/Uni/Eagle_Master/Hiwijob_IZW/Remoteness_tutorial"
# For GADM download set name, iso code and level for required country
name <- "cambodia"
iso <- "KHM"
level <- 1
# Set buffer that should be applied to aoi
buffer <- 10000
# Set EPSG code of a metric CRS (e.g utm)
crs_m <- 32648
Download country data.
GADM provides maps and spatial data for all countries and their
sub-divisions. You can download your required country data from inside
R. Alternatively download from website.
# Get or import country data
country <- geodata::gadm(iso, level=level, path = tempdir())
country <- st_as_sf(country)
# # Export to folder
# st_write(country,
# dsn = file.path(wd, "data/gadm/"),
# layer = paste0("gadm_", name, "_", level),
# driver = "ESRI Shapefile")
Select or import aoi.
In this tutorial the province of Kâmpóng Thum is chosen as aoi.
Define/import your aoi for which you want to calculate remoteness.
# Select province as example aoi, or import your Area of Interest
province <- "Kâmpóng Thum"
aoi <- country[(country$NAME_1 %in% province),]
# Export to folder
# st_write(aoi,
# dsn = file.path(wd, "data/gadm/"),
# layer = "aoi",
# driver = "ESRI Shapefile")
mapview(country, layer.name = name) + mapview(aoi, col.regions = "orange")
Figure: Interactive map of study area in cambodia.
Download OpenStreetMap (osm) data by country.
Download data in osm.pbf format!
PBF files are generally much smaller than OSM XML or SHP files.
OpenStreetMap Protocolbuffer Binary Format (PBF) is an open source transfer format for vector GIS data created by the OpenStreetMap community. More info here. A provider of osm data is the Geofabrik Download Server.
Options 1&2: Either you download the osm data on a country-based level and optionally clip it to your buffered study area in R …
Option 3: … or you download the data with the QGIS plugin QuickOSM (only possible for small study areas, needs to be tried out).
Option 1: download osm.pbf data from Geofabrik with
osmextract package in R.
# Insert name of country or insert aoi (sf file!)
its_details <- oe_match(name) # alternative: oe_match(aoi)
## The input place was matched with: Cambodia
# Download data to folder
its_pbf <- oe_download(
file_url = its_details$url,
file_size = its_details$file_size,
#provider = "test",
download_directory = file.path(wd, "data/osm")
) # --> Time difference of 7.108055 secs
## The chosen file was already detected in the download directory. Skip downloading.
Option 2: download osm.pbf data from Geofabrik website and import to R.
Option 3: download osm.pbf data via QuickOSM plugin in QGIS.
Import duration depends on size. For cambodia it takes ~ 30 seconds, for larger data sets it can take up to several minutes.
# list files in folder
list.files(file.path(wd, "data/osm"))
# Import osm.pbf file
roads <- st_read(# filename downloaded with osmextract (option 1)
dsn = file.path(wd, "data/osm/geofabrik_cambodia-latest.osm.pbf"),
# filename downloaded from website (option 2)
#dsn = file.path(wd, "data/osm/cambodia-latest.osm.pbf"),
query = "select highway from lines")
Select only drivable road categories. More info on specific road categories at https://wiki.openstreetmap.org/wiki/Key:highway
Check for NA.
# Check for NA
any(is.na(roads))
## [1] TRUE
# If NA is present, rename to "unknown" and optional check whether to include or not
roads$highway[is.na(roads$highway)] <- "unknown"
# # Optional map or export only "unknown" roads to check if they are drivable or not
# na <- roads[roads$highway %in% "unknown",]
# mapview(na) # -> check in R
# st_write(na, file.path(wd, "data/osm/na_unknown_roads.shp")) # -> check e.g in QGIS
Select drivable road categories.
# List single road categories in alphabetical order
sort(unique(roads$highway))
## [1] "bridleway" "construction" "corridor" "cycleway"
## [5] "disused" "elevator" "footway" "living_street"
## [9] "motorway" "motorway_link" "path" "pedestrian"
## [13] "primary" "primary_link" "proposed" "raceway"
## [17] "residential" "rest_area" "road" "secondary"
## [21] "secondary_link" "service" "services" "steps"
## [25] "tertiary" "tertiary_link" "track" "trunk"
## [29] "trunk_link" "unclassified" "unknown" "unused"
# Select only drivable roads (check if is complete, "unknown" roads were NOT included in this tutorial)
include <- c("living_street", "motorway", "motorway_link", "primary", "primary_link", "proposed", "residential", "road", "secondary", "secondary_link", "tertiary", "tertiary_link", "trunk", "trunk_link", "unclassified")
roads <- roads[roads$highway %in% include,]
First transform to metric coordinate reference system.
Then apply buffer to aoi and clip roads to buffered aoi.
# Convert to metric CRS
country_m <- st_transform(country, crs_m)
aoi_m <- st_transform(aoi, crs_m)
roads_m <- st_transform(roads, crs_m)
# Clip roads to buffered aoi
roads_m_aoi <- st_intersection(roads_m, st_buffer(aoi_m, buffer))
# Export to folder
# st_write(roads_m_aoi,
# dsn = file.path(wd, "data/osm/"),
# layer = "aoi_buff_drivable_roads",
# driver = "ESRI Shapefile")
# Map roads
mapview(country_m, layer.name = name, col.regions = "black", alpha.regions = 0) +
mapview(aoi_m, layer.name = "aoi", col.regions = "orange") +
mapview(roads_m_aoi, layer.name = "osm drivable roads")
info: if other starting points …
Osm roads are multiline features, but the algorithm that calculates remoteness requires starting points. Therefore, lines need to be converted to equal spaces points. Using the “Points along geometry” tool in QGIS (v.3.16.10) is the easiest and fastest way to perform this step. The distance between the points was set to 100 meter in this tutorial.
Google Earth Engine (GEE) is a cloud computing platform, powered by Google cloud infrastructure. It offers new opportunities for Big Data analyses with huge data sets up to petabyte scale, and analyses are run on Google servers, which have much higher computational power and storage capacity than most local computer sources.
The remoteness analysis can only be performed in GEE. Users must log in with a Google account in order to work with the user interface (uploading input file, running the script, exporting products). Hence it is not open source software but costfree.
More information on how to create a Google Earth Engine account can
be found here.
To get familiar with the user interface check out this
page.
Useful information on script and analysis performance in GEE:
Follow this link to open the remoteness analysis script in GEE.
Open script link and “Run” demo to see how script performs remoteness analysis with example aoi and access points.
GEE script variables: geometry, startpoints
Before working with own data, comment out demo parameters!
GEE script variables: geometry, startpoints
Upload area of interest and access points (previously created in R / QGIS) as shapefiles in ‘Assets’ tab.
Uploaded files will appear soon in ‘Assets’ tab (left window, dotted magenta line). If not click “refresh” (solid magenta line). Copy the asset paths into script (dotted magenta lines in code editor).
Optionally draw aoi as polygon on map. The file will automatically appear as “geometry” variable on top of the script. In this case make sure all other “geometry” variables in the script (demo, asset import) are commented out.
GEE script variables: watermask, occ, buffer, maxPixels
Optionally change the buffer to be applied around aoi, but it should be the same buffer that was already applied in previous steps when clipping roads.
Only modify maximum pixels if export fails. See section 4.5.
Finally click “Run” at the top of the script (solid magenta line).
After executing the script, some of the selected parameters will be displayed in the print console (right window, dotted magenta line).
Study area, starting points, cost raster and the cumulative cost raster (remoteness layer) are added to the map (hover over layer panel box in map). Optionally modify visualization parameters of single layers When clicking on “settings” of respective layer, or uncheck layer if it should not be displayed in map.
Refine water mask
If a water mask is applied (var watermask = “T”), 10 different cost raster are added to the map by default, based on different values for the probability of water occurrence (10%, 20%, … 100%).
This way the layers can be compared with each other (check / uncheck). If the water mask should be refined, the occ parameter must be modified and the script must be re-executed. The water mask applied to the images ready for export are based on the specified occurrence value.
Figure: Compare the output rasters of the remoteness analysis after applying different water masks.
Final image export needs to be initiated in ‘Tasks Tab’ (right
window).
1. Go to “Tasks” tab.
-> Optionally modify default export parameters (dotted magenta line).
Depending of size of study area, an image export can take from minutes to days. Multiple exports can be run in parallel.
For example, if maxPixels parameter would have been set to 4000, an error message would occur after initiating the final export task (see Figure). In this case, the maxPixels value (see section 4.3) would have to be increased (see error message) and the script and export task would have to be executed again.